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I.  INTRODUCTION 


A  current  focus  of  interest  in  the  U.8.  Any  ballistic  research  program 
involves  the  numerical  calculation  of  compressible  flow  about  muz ale  brake 
devices.  By  absorbing  a  portion  of  the  recoil  impulse,  the  muzzle  brake 
penits  design  of  large  caliber  weapons  characterized  by  increased  range  with¬ 
out  increased  weight.  Near  field  calculations  are  helpful  in  studying  fatigue 
life  and  structural  integrity  of  blast-loaded  surfaces.  Far-field  calculations, 
which  often  can  be  performed  with  less  complex  flow  models,  can  detenine  if 
safe  maximum  overpressures  exist  in  the  gun  crew  area. 

Near-to-intenediate  range  muzzle  brake  calculations  have  recently  been 
obtained1  using  a  popular,  locally  one-dimensional,  first-order  accurate  method 
of  Godunov.2  Model  results  afford  predictions  of  peak  over-pressure  levels  in 
the  proximity  of  the  brake  and  provide  initial  data  for  continuing  far-field 
calculations  by  independent  means. 

In  the  interest  of  economizing  computer  resources,  for  far-field  calcula¬ 
tions,  it  is  desirable  to  employ  simple  one-dimensional  flow  models.  Mach 
contours  for  the  early  stages  of  a  typical  blast  are  exhibited  in  Figure  1. 
Already,  local  spherical  symmetry  is  suggested,  and  evidence  from  spark 
photography  confirms  that  the  trend  is  characteristic  of  later  evolution.1 
Under  this  hypothesis,  some  exploratory  far-field  calculations  are  underway, 
which  utilize  spherically  symmetric  flow  models  and  the  numerical  method  of 
characteristics.  3»H 


1.  G.  E.  Widhopf ,  J.  C.  Buell,  and  E.  M.  Schmidt ,  ' Time  Dependent  Bear-Field 
Muzzle  Brake  Simulations,"  AIAA-82-0973,  AIAA/ASME  3rd  Joint  Thermo- 
physics  ,  Fluids  and  Heat  Transfer  Conference,  St.  Louis,  Missouri,  June 
1982. 

2.  A.  M.  Godunov,  A.  V.  Zabrodin,  and  G.  P.  Prokopov ,  "Difference  Schemes  for 
Two-Dimensional,  Unsteady  Problems  in  Gas  Dynamics  and  Calculation  of  Flows 
With  a  Detached  Shock  Wave, "  Journal  of  Computing  Mathematics  and  Mathe¬ 
matical  Physics,  USSR  Academy  of  Sciences,  Vot.  1,  No.  6,  November  - 
December  1961.  (Translation) 

3.  M.  L.  Bundy,  "A  Nonsimilar  Solution  for  Blast  Waves  Driven  by  an  Asymp¬ 
totic  Piston  Expansion, "  AIAA-83-0496 ,  AIAA  21st  Aerospace  Sciences 
Meeting,  January  1983,  Reno,  Nevada  and  U.  S.  Army  Ballistic  Research 
Laboratory ,  Aberdeen  Proving  Ground,  MD,  BRL  Technical  Report  ARBRL-TR- 
02497,  June  1983.  (AD  A130012) 

4.  M.  L.  Bundy,  C.  H.  Cooke,  and  E.  M.  Schmidt,  "Reshaping  an  Artificially 
Diffused  Shock  Solution, "  BRL  Report,  to  be  published. 


The  objective  of  the  present  research  is  to  investigate  an  alternative 
numerical  method ,  which  might  complement*  or  perhaps  supplant*  the  method  of 
characteristics  approach.  Our  intent  is  to  revise  Van  Leer's  second-order 
accurate*  Godunov  method, 5  originally  formulated  in  the  framework  of  Lagrangean 
fluid  dynamics*  for  application  to  the  one-dimensional  Euler  equations.  A 
similar  effort,  as  yet  unpublished*  and  which  differs  somewhat  in  philosophy, 
has  been  carried  out  by  Colella. 

The  desirability  of  investigating  second-order  accurate  methods  for  shock 
capturing  is  illustrated  by  Figure  2.  Here*  the  first-order  Godunov  method* 
implemented  by  personnel  of  Aerospace  Corporation  for  numerical  solution  of 
the  two-dimensional  axis-symmetric  Euler  equations*  has  been  applied  to  calcu¬ 
late  a  conical  flow  which  simulates  some  of  the  more  dominant  characteristics 
of  a  muzzle  blast.  Assuming  there  is  provided  some  heuristic  model  of  contact 
surface  history*  as  well  as  initial  data  between  the  outer  shock  and  its  driving 
contact  surface*  this  calculation  could  be  continued  into  the  far-field  by  the 
method  of  characteristics. 3  However,  shock  smearing  due  to  the  artificial  vis¬ 
cosity  of  the  numerical  method  in  this  case  makes  troublesome  the  question  of 
precise  shock  location  and  strength.  A  reshaping  of  initial  data  near  the  shock* 
or  else  a  more  accurate  calculation  which  provides  crisper  shock  structure* 
appears  to  be  called  for. 

In  the  past  few  years,  a  variety  of  new  methods  for  numerical  calculation 
of  flows  with  embedded  shocks  has  evolved,  of  which  references  5-9  are  perhaps 
a  representative  sample.  Van  Leer's  second-order  sequel  to  the  original  Godunov 
siethod  appears  among  the  more  promising.  The  method  is  alledged5to  give 


5.  B.  Van  Leer ,  "Towards  The  Ultimate  Conservative  Difference  Scheme.  V. 

A  Second-Order  Sequel  to  Godunov  ’e  Method , "  Journal  of  Computational 
Physics,  Vol.  32,  pp.  101-136 ,  1979. 

6.  P.  Colella,  "A  Direct  Eulerian  MUSCL  Scheme  for  Gas  Dynamics, "  Lawrence 
Berkeley  Laboratory  Report  LBL-14104,  February  1982. 

7.  J.  L.  Steger  and  W.  F.  Warming ,  "Flux  Vector  Splitting  of  the  Inviscid 
Gas  Dynamics  Equations  with  Application  To  Finite  Difference  Methods, " 
Journal  of  Computational  Physics,  Vol.  40,  No.  2,  April  1981. 

8.  G.  Moretti,  "The  \  Scheme , "  Computers  and  Fluids.  Vol.  7,  pp.  191-205, 
Pergamon  Press,  1979. 

9.  H.  C.  Tee  and  R.  J.  Warming,  "Implicit  Total  Variation  Diminishing 
Schemes  For  Steady  Flow  Calculations ,"  AIAA-83-1902,  AIAA  6th  Compu¬ 
tational  Fluid  Dynamics  Conference,  Danvers,  Massachusetts ,  July 
1983. 
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superior  resolution  of  shocks  and  flow  discontinuities,  compared,  say,  to  the 
methods  surveyed  by  Sod10  and  Miner  and  Skop.11  However,  remapping  from  the 
Lagrangean  to  Euler  variables  requires,  perhaps  significant,  extended  comput¬ 
ing  time  per  cycle.5  Hence,  it  appears  an  unnecessary  encumbrance. 

The  purpose,  then,  of  this  research  is  to  reformulate  Tan  Leer's  algo¬ 
rithm  in  Eulerian  fluid  dynamics  framework,  revising  as  it  becomes  necessary, 
in  order  to  achieve  a  more  accurate,  one-dimensional  shock  capturing  algo¬ 
rithm,  which  could  also  be  employed  in  two-dimensional  calculations  through 
fractional  splitting  of  the  equations  of  flow. 

II.  GOVERNING  EQUATIONS 

In  strong  conservation  law  form,  the  Euler  equations  for  one-dimensional 
ideal  compressible  flow  can  be  written 


F  *  P  +  pu 
(P  +  pE)u 


(P+pE)u 


where  o  "0,1,  or  2,  in  turn  for  cartesian,  cylindrical,  or  spherical  coordi¬ 
nates.  R  is  the  respective  distance  coordinate. 

The  fluid  dynamic  variables  are: 

c  ■  local  speed  of  sound  in  fluid 
p  ■  density 
u  •  velocity 

10.  G.  A.  Sod ,  " A  Survey  of  Several  Finite  Difference  Methods  For  Systems  of 
Hyperbolic  Conservation  Lots.”  Journal  of  Computational  Physics ,  Vol.  27, 
pp.  1-31 ,  1978. 

11.  E.  W.  Miner  and  S.  A.  skop ,  ” Explicit  Time  Integration  For  The  Finite 
Element  Shock  Wave  Equations,"  AIAA-8 2-0994,  ASME/AIAA  3rd  Joint  Thermo¬ 
physics,  Fluids,  Plasma  and  Heat  Transfer  Conference,  St.  Louis,  Missouri, 
June  1982. 
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p  “  pressure 

E  *■  specific  total  energy 
e  ■  specific  internal  energy. 
Here  Cp, cy  are  specific  heats,  and 

Y  -cp/cv 
c2  -  yP/  p 

P  ■  ( y  ~  1)  Pe 


E  ■  e  ♦  u 

2 


III.  GODUNOV  METHODS 

We  shall  derive  the  Godunov  algorithm  for  the  case  of  a  uniform  grid; 
however,  the  method  is  readily  adaptable  to  encompass  non-uniformity.  By 
integrating  Equation  (1)  over  a  typical  space-time  cell  h  * R  f R  itiL 
tjj  <  t<tn+i  and  applying  Green  s  theorem  for  the  plane,  ve  arrive  at  the 
exact  equation 


i  - 1 


ui.i 


1  *  1  ,‘n  .  1 

M<f>  -  —  I 

AR  AR  I 

i  ^  t 


R.  . 

*1  +  1 


G  dR  dt. 


Here,  superscript  usage  of  a  space  index  denotes  advanced  time  level, 
while  corresponding  subscript  usage  denotes  present  time  level:  i.e., 

f  '  <  h,n  *  1 

(  >1  '  <  >i.n. 

The  space  average  over  a  cell,  space-centered  at  R  ^  +  j  is 


Ri  +  1 


“i  *  i  ■  m  I 


U(R,t  )  dR,  AR  *  R.  .  -  R.; 
n  l  +  l  i 


'I. 


*1 


while  the  time  average  flux  on  interface  R  ■  Rj  ii 


F(U(R.,t))dt,  At  «  tn  +  i  -  tn.  (6) 
t 

n 


Godunov's  method  can  be  made  first- or  second-order  accurate*  depending  upon 
hov  the  flux  integrals*  Equation  (6)  and  the  cell  integral  of  G  in  Equation 
(3)  are  approximated. 


A.  First-Order  Accurate  Hethod 


In  Godunov's  original  derivation,2  the  function  U(R  ,tn)  is  approxi¬ 
mated  with  a  piecewise  constant  function  which  on  cell  Rj  <.  R  <  R^  +  j 
assumes  the  average  value  0^ +  j.  Cell  averages  are  updated  to  the  next  time 
level  through  approximating  the  integrals  in  Equation  (3),  by  solving  a 
Riemann  problem  at  each  cell  interface.  The  Riemann  problem  entails  the 
solution  of  Equation  (1)  for  t>t  ,  with  am  0  and  with  initial  conditions  at 
t  •  t  given  by:  n 


Vi  ♦  l. 


'i  >  i. 


R  >R± 

R  <R. 
x 


(7) 


where  V  •  (P,  p,u).  However*  only  the  resulting  values  of  V  on  the  inter¬ 
face  Rj.  are  of  interest.  The  primary  mechanism  for  interaction  of  the  dis¬ 
continuity,  Equation  (7),  i*  the  propagation  of  expansion  or  compression  waves 
away  from  the  interface.  The  nonlinear  equations  which  afford  an  iterative 
solution  of  the  Riemann  problem  are  well-documented  in  references  5,  10,  and 
12. 


For  purposes  of  numerical  stability*  the  time  step  At  is  chosen  to  be 
such  that  propagation  times  are  insufficient  to  allow  waves  from  adjacent 
discontinuities  to  have  influence  on  interface  R  ■  Rj..  Then  values 


*  lim  V(Ri,t) 


t  -*•  t 


n  * 


(8) 


12.  M.  Holt t  numerical  Fluid  Dunamice,  Springer-Verlagt  Berlin ,  Heidelburg , 
New  York ,  1977. 


obtained  from  solving  the  Riemann  problem  at  each  interface,  together  with 
cell  average  values,  are  used  to  approximate  the  integrals  in  Equation  (3). 
Godunov's  first-order  accurate  method  results: 


_  i  +  1 

U 


U. 


At 


i*rs  lF(ui  ♦  i> 


-  F(U*)]  -  At  G(U.  +  j). 


(9) 


B.  Second  Order  Method 


For  the  second-order  accurate  Godunov  method,  primitive  quantities 
stored  at  each  time  level  are  cell  average  Ui  +  \  and  interface  differences 
IV]^  +  |  ■  Vi  +  i  -  V^.  The  calculation  of  the  average  derivative  is  now 
possible: 


i  ♦  l 


9 V 

9Ri  +  1 


1_ 

AR 


[V] 


i  ♦  1 


AR 


(10) 


This  affords  a  more  accurate,  piecewise  linear  function  approximation:  On 
a  cell  R.  <  R  <  R.  +  lt 


V  =  Vi  ♦  1  +  (Vi  ♦  1 


(R  -  R, 


i  +  lJ 


(11) 


Corresponding  inputs  for  the  Riemann  problem  at  interface  R^  are,  from 
Equation  (11), 


Vi  +  ■  Vi  ♦  !  *  iIV>i  .  1  •  tl2> 

The  output  from  the  Riemann  problem,  solved  as  previously,  is  the  value 


V.  =  lim  V(R.,t), 
t  -►  t* 


(13) 


which  results  immediately  after  the  interface  discontinuity  is  resolved.  In 
addition,  a  value 


<  -  ft  <Ri-« 


(14) 


for  the  corresponding  resolved  time  derivative  is  to  be  obtained  by  auxiliary 
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means  (see  next  section).  As  established  by  Van  Leer*1  the  interface  approxi' 
nation 


v  =  V*  ♦  (VJ  *  (t  -  t J  +  0  (At)2 


(15 


can  be  applied  to  evaluate,  with  a  higher  order  of  accuracy,  the  flux 
integrals  in  Equation  (3).  The  cell  integral  can  be  approximated,  employing 
the  trapezoidal  rule,  with 


G  dR  dt  = 


(16 


[G(R.  +  j,t)  +  G(R.,t)]dt  -  0  (AR)3 


which  now  involves  interface  values.  For  a  typical  interface  integral 


1 

G(U,R.)dt 


[G(U\R.)  +  GfUj.R^]  +  0(At)3. 


Advanced  time  level  interface  values  are  predicted  by  means  of 


V1  =  V*  +  (V,)*  At  +  0  (At)2, 

l  t  l 


(17 


(18 


and  initial  interface  differences  with 

[V]i  *  *  -  V1  *  1  -  V*. 


(19 


At  the  cost  of  additional  storage  and  altered  processing  stream,  after 
solution  of  the  Riemann  problem  adjusted  interface  differences 


m 


could  be  computed,  to  be  used  for  more  accurate  evaluation  of  the  interface 
time  derivatives,  described  below.  However,  Van  Leer  does  not  appear  to 
have  used  this  device,  and  we  have  not  verified  whether  the  additional  cost 
is  justified. 

C.  Resolved  Interface  Initial  Time  Derivatives 

For  the  compressible  flow  equations.  Equations  (1-2),  compatibility 
relations  along  characteristics  are. 


dP  .  2  dp 
dt  =  c  Ht 


dR 

0n  dt  =  U 


du  +  1_  dP  _  +  cruc 
dt  ”  pc  dt  R 

dR 

on  =  u  ±  c. 


(21) 


(22) 


To  better  insure  correct  transmission  of  signals,  these  Equations  are 
differenced  (spatially)  as  in  Moretti's  A  -  scheme  ,B  in  order  to  obtain 
relations  which  can  be  solved  for  initial  interface  time  derivatives.  This 
differencing  is  given  by  the  equations 


*  *  ★ 

—  +  (i— )  (— ■)  s 

3t.  V.  l3tJ. 


,  ..  r3u  1  9P,  A  ,ouc. 

•  (u  *  c)i*  hr’  j 


(23) 
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r-V 
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- 

i  i  i 


,  .  r9u  1  3P-.  ,ouc.. 

-  tu-c)  i-{3R-pE3R}  .  -  hr >. 

1-  l 


(24) 


l£.  =  1 _  2  {—  +  (u  — )  }  -  (u  — ) 

9t .  (c*.r  l9t.  lU  1  9R; 


(25) 


it 


14 


it 


Where  space  derivative  occur  in 
are  used,  on  the  (±)  aide  of  an 
u  is  positive,  up-  or  dovn-vind 


Equations  (23-29),  average  space  derivatives 
interface.  Depending  upon  whether  or  not 
differencing  is  employed  in  Equation  (25). 


IV.  COMPUTATIONAL  RESULTS 

The  first  and  second-order  Godunov  methods  previously  discussed  have 
been  applied  to  the  linear  shock  tube  calculation  reported  by  Miner  and 
Shop.11  Here,  an  infinite  tube  contains  gas  in  two  compartments  initially 
separated  by  a  diaphragm.  Table  I  shows  the  respective  initial  conditions. 


TABLE  1.  LINEAR  SHOCK  TUBE  DATA 


P 

o 


P 


o 


1. 

1. 


u 


o 


-  0. 


.1 


Pj  -  .125 


u. 


-  0. 


The  Godunov  calculation  is  programmed  to  choose  its  own  time  step,  in 
accordance  with  stepwise  stability  restrictions.  The  results  after  one 
hundred  cycles,  for  the  first-order  accurate  calculation,  are  displayed  in 
Figures  3-6.  Figures  4-5,  in  particular,  show  the  smearing  of  shock  struc¬ 
ture  due  to  the  inherent  numerical  dissipation  of  the  method,  present  even 
on  this  very  fine  grid. 

The  second-order  accurate  method  is  activated  by  a  program  switch. 
Figures  7-10  show  the  results  after  another  one  hundred  cycles  of  calcula¬ 
tion.  An  immediate  sharpening  of  the  shock  is  to  be  observed. 


It  seems  to  be  a  consensus  of  opinion  that  higher  order  methods  for 
shock  capturing  are  likely  to  be  characterised  by  overshoot  and  oscillations 
behind  the  shock.13  Our  results  appear  to  be  no  exception.  Van  Leer's 
oscillation  limiting  techniques5  were  attempted,  as  well  as  sparse  use  of 
numerical  viscosity  in  the  vicinity  of  the  shock.  For  our  code,  the  second 
approach  seemed  to  give  as  good  results  as  the  first.  Here,  the  Riemann 
solver  provides  a  shock  Mach  number,  which  is  a  maximum  at  the  point  of  in¬ 
flection  occurring  within  the  structure  of  the  physical  shock.  This  pro¬ 
vided  a  means  for  limiting  application  of  artificial  viscosity,  to  a  few 
points  either  side,  but  concentrated  more  to  the  upstream  side  of  the  shock. 
For  density  and  velocity,  the  artificial  viscosity  was  chosen  as 


V  *  10 


3 

Pi  ♦  1 

-  2P. 

l 

*Pi  -  ! 

Pi  ♦  1 

+  2P. 

l 

+  Pi  -  1 

(vi .  r2  vi  ♦  vi  - 1);  (26) 


13.  J.  L.  Sieger ,  private  coirrmnioation. 
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The  result*  of  Figures  7-10  are  calculated  vith  this  dissipation,  which 
is  of  the  order 


V  =  0  (10-6  32V) 


The  effects  seem  to  be  a  dampening  of  oscillations  behind  the  shock,  with 
no  apparent  degradation  in  crispness.  However,  the  overshoot  appears  to 
persist,  about  3Z  in  error. 

In  order  to  see  that  the  shock  is  propagating  properly,  the  calculation 
with  the  second-order  method  has  been  allowed  to  continue  to  time  t  -  .14.  As 
reported  by  Miner  and  Skop, 11  at  this  time  the  shock  front  should  have  pro¬ 
gressed  to  R  ■  .75.  Figures  11-17.  verify  that  the  shock  front  is  approxi¬ 
mately  at  this  location. 


V.  SUMMARY  AND  CONCLUSIONS 

First  and  second-order  accurate  Godunov  methods  for  the  numerical  solu¬ 
tion  of  ideal  compressible  flows  with  embedded  shock  waves  have  been  formu¬ 
lated,  programmed,  and  tested  by  means  of  a  linear  shock  tube  calculation. 
Results  show  an  immediate  improvement  of  the  crispness  in  shock  structure, 
for  the  higher  order  accurate  method.  The  higher  order  method  appears  to  have 
inherent  overshoot  at  the  shock,  together  with  oscillations  behind  the  shock. 
It  appears  mandatory  to  dampen  these  oscillations,  by  means  of  Van  Leer-type 
oscillation  limiting  schemes,  or  else  through  addition  of  artificial  viscos¬ 
ity  restricted  to  a  small  region  behind  the  shock.  Surprisingly  enough,  for 
the  present  problem  both  schemes  were  found  to  give  comparable  results  in 
this  respect. 

Perhaps  we  should  mention  how  the  present  method  differs  from  Van  Leer's 
original  version, 5  aside  from  the  conversion  to  Eulerian  fluid  dynamics. 

Major  differences  are: 

a.  Omission  of  some  near-shock  terms  from  Equations  (23-24),  heuristi- 
cally  added,  perhaps  to  insure  entropy  increase  at  the  shock. 

b.  Use  of  Godunov's  original  nonlinear  iteration  scheme  (References  2, 
12)  for  the  Riemann  problem,  versus  Van  Leer's  more  elaborate  accelerated 
version.  Although  slower  computationally,  this  scheme  does  distinguish  be¬ 
tween  shock  and  compression  waves;  hence, it  should  be  more  accurate,  as 
evidenced  by  the  omission  of  (a). 

c.  For  the  Lagrangean  formulation,  the  contact  surface  in  the  Riemann 
problem  lies  along  the  cell  interface.  Hence,  it  is  reasonably  for  Van  Leer 
to  employ  separate  left  and  right  density  Pj+  in  calculating  Pg  .  However, 
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for  the  Euler  version  this  practice  does  not  appear  as  logical,  since  the 
contact  surface  can  lie  on  either  side  of  the  cell  interface.  The  reault, 
and  possibly  the  biggest  drawback  of  the  Euler  version,  ia  that  contact  surface 
resolution  does  not  such  improve,  when  going  over  to  second-order  accuracy. 


Smearing  Typical  of  a  First  Order  Accurate  Godunov  Calculation 


:irst  Order  Godunov  Calculation 


First  Order  Godunov  Calculation 


Figure  6.  First  Order  Godunov  Calculation 
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Second  Order  Godunov  Calculation 


Progression  of  Pressure  Jump,  Second  Order  Godunov 
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